function lli= indiv_likelihood(lambda,  rho, Deltas, default, probmf,c1, c2)
% this function returns a vector with the individual log likelihoods

% inputs: 
% lambda - parameters to be estimated
% rho - vector with observations of donations
% Deltas - Delta: vector with utility gain of donating rho instead of the
% default
% default - the default value the agent had (10, 20 or 50)
% probmf - relative frequency of bins of donations from 0 to 300
% c1 - rescaling parameter of lambda2: lambda2/c1
% c2 - rescaling parameter of lambda3: lambda3/c2


%-------unpack STATS
lambda1= lambda(1);
lambda2a= lambda(2);
lambda2=lambda2a/c1;
lambda3a=lambda(3);
lambda3=lambda3a/c2;
N=length(rho);
% preallocate variables
lli=zeros(N,1);

% assign default values
d1=10; d2=20; d3=50;


% assign probability mass to each observed value of rho
f=estimpmf(probmf, rho);

% calculate sums for x=0 and x=d in probability (functions defined at the bottom)
Sum0_d1=Sum_lik0(lambda2, lambda3, d1, probmf);
Sum0_d2=Sum_lik0(lambda2, lambda3, d2, probmf);
Sum0_d3=Sum_lik0(lambda2, lambda3, d3, probmf);

Sumd_d1= Sum_likd(lambda2, lambda3, d1, probmf);
Sumd_d2= Sum_likd(lambda2, lambda3, d2, probmf);
Sumd_d3= Sum_likd(lambda2, lambda3, d3, probmf);

%-- log likelihood for each observation i up to N; define cases as in
%Appendix section D

%- case 1: x=0
case1_d1= (default==d1 & rho==0);
case1_d2= (default==d2 & rho==0);
case1_d3= (default==d3 & rho==0);
%- case 2: x=d
case2_d1= (default==d1 & rho==d1);
case2_d2= (default==d2 & rho==d2);
case2_d3= (default==d3 & rho==d3);
%- case 3: 0<x<d/2
case3_d1= (default==d1 & rho<(d1/2) & rho>0);
case3_d2= (default==d2 & rho<(d2/2) & rho>0);
case3_d3= (default==d3 & rho<(d3/2) & rho>0);
%- case 4: x>=d/2, x~=d
case4_d1= (default==d1 & rho>=(d1/2) & rho~=d1);
case4_d2= (default==d2 & rho>=(d2/2) & rho~=d2);
case4_d3= (default==d3 & rho>=(d3/2) & rho~=d3);

%- case 1 individual log likelihood
lli(case1_d1)=log(f(case1_d1) + (1-lambda1)*Sum0_d1);
lli(case1_d2)=log(f(case1_d2) + (1-lambda1)*Sum0_d2);
lli(case1_d3)=log(f(case1_d3) + (1-lambda1)*Sum0_d3);

%- case 2 individual log likelihood
lli(case2_d1)= log(f(case2_d1) + (1-lambda1)*Sumd_d1);
lli(case2_d2)=log(f(case2_d2) + (1-lambda1)*Sumd_d2);
lli(case2_d3)=log(f(case2_d3) + (1-lambda1)*Sumd_d3);

%- case 3 individual log likelihood
% let z be the second term of the probability, 
% function case3_lik is defined at the bottom
z=case3_lik(lambda1, lambda2, lambda3, default, Deltas, rho);

lli(case3_d1)= log(f(case3_d1)) + log(z(case3_d1));
lli(case3_d2)=log(f(case3_d2)) + log(z(case3_d2));
lli(case3_d3)=log(f(case3_d3)) + log(z(case3_d3));

%- case 4 individual log likelihood
lli(case4_d1)= log(f(case4_d1)) + log(lambda1 + (1-lambda1).*(1-exp(-lambda2.*Deltas(case4_d1))));
lli(case4_d2)=log(f(case4_d2)) + log(lambda1 + (1-lambda1).*(1-exp(-lambda2.*Deltas(case4_d2))));
lli(case4_d3)=log(f(case4_d3)) + log(lambda1 + (1-lambda1).*(1-exp(-lambda2.*Deltas(case4_d3))));

end

function Sum0= Sum_lik0(lambda2, lambda3, d, probmf)
% this function calculates the sum in the probability for case 1
% sum over rho, 0<rho<d/2
r=(1:(d/2 -1))' ;
def= d.*ones(length(r),1);
% probability mass function: elements correspond to rho-1 
pmass=probmf(2:(d/2));

Del= ((r.^2)+(def.^2))./2 - r.*(def); 
Sum0_i=  pmass.*((lambda3/(lambda2+lambda3)).* exp(-lambda2*(r.^2)./2)-...
    (lambda3/(lambda2+lambda3)).*exp(-lambda3*((def.^2)/2 - r.*def)).*exp(-lambda2.*Del)); 
Sum0=sum(Sum0_i); 

end

function Sumd= Sum_likd(lambda2, lambda3, d, probmf)
% this function calculates the sum in the probability for case 2
% rho from 1 to 300
r=(1:300)';
def= d.*ones(length(r),1);
% probability mass function: elements correspond to rho-1 
pmass=probmf(2:301);
Sumd_i=zeros(length(r),1);
Del= ((r.^2)+(def.^2))./2 - r.*(def); 
% sum for rho <= d/2
a=(r<=(d/2));
% sum for rho > d/2 & rho~=d
b=(r>(d/2) & r~=d);
 
Sumd_i(a)=  pmass(a).*exp(-lambda3.*((def(a).^2)./2 - r(a).*def(a))).*exp(-lambda2.*Del(a));
Sumd_i(b)= pmass(b).*exp(-lambda2.*Del(b));
Sumd=sum(Sumd_i);
end

function z=case3_lik(lambda1, lambda2, lambda3, default, Deltas, rho)
% this function calculates the second term of the probability for case 3
% assign default values
d1=10; d2=20; d3=50;
% define case 3
case3=((default==d1 & rho<(d1/2) & rho>0) | (default==d2 & rho<(d2/2) & rho>0)| (default==d3 & rho<(d3/2) & rho>0));
% preallocate
z=zeros(length(rho),1);
a=zeros(length(rho),1);
a(case3)=-lambda2.*Deltas(case3) - lambda3.*(((default(case3)).^2)./2-rho(case3).*default(case3));
z(case3)= lambda1 + (1-lambda1).*(1-exp(-lambda2.*((rho(case3)).^2)./2)  + (lambda2/(lambda2+lambda3)).*(exp(-lambda2.*((rho(case3)).^2)./2) - exp(a(case3))));
end

